Universality in the vibrational spectra of weakly-disordered two-dimensional clusters 
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We numerically investigate the vibrational spectra ol single-component clusters in two-dimensions. 
Stable configurations of clusters at local energy minima are obtained, and for each the hessian matrix 
is evaluated and diagonalized to obtain eigenvalues as well as eigenvectors. We study the density of 
states so obtained as a function of the width of the potential well describing the two-body interaction. 
As the width is reduced, as in three dimensions, we find that the density of states approaches a 
common form, but the two-peak behavior survives. Further, calculations of the participation ratio 
show that most states are extended, although a smaller fraction of the degrees of freedom are involved 
in these modes, compared to three dimensions. We show that the fluctuation properties of these 
modes converges to those of the Gaussian orthogonal ensemble of random matrices, in common with 
previous results on three dimensional amorphous clusters and molecular liquids. 



PACS numbers: 63.50.-x, 24.60.Lz, 05.45.Mt, 61.43.Fs 



INTRODUCTION 

The study of vibrational spectra has long provided im- 
portant insights into the understanding of amorphous 
states of matter [H, such as molecular glasses. However 
many key questions still remain unanswered. The vibra- 
tional disorder that is present in systems such as glasses 
is mainly topological in origin and over the years, a lot 
of effort has been put in to study glasses and glass tran- 
sitions in terms of underlying potential energy surface 
(PES) 0, S, i, i, S, 0, &, la]. One of the main focus has 
been the study of inherent structures, i.c local minima 
of PES and properties like aging have also been stud- 
ied as a function of inherent structures [l3|- Regarding 
vibrational spectra of systems having a topological disor- 
der 0,|il[ii,[ii,[il,|iali3,S[2i|,El[2i,[2i,M 



28[, a popular approach rely on the study of statistical 
properties of random matrices, as within the harmonic 
approximation, the dynamical or the Hessian matrix con- 
tains all the dynamical features 11, 13, 13, 3, 15, 3, 13] ■ 



Recently, studies of vibrational spectra of amorphous 
systems have revealed some interesting universal aspects. 
In particular, studies of three-dimensional amorphous 
clusters [l^ 0, Ell have demonstrated the universality 
in the functional form of the density of states (DOS) 
and further, the statistical fluctuations have been shown 
to obey the characteristics of the Gaussian Orthogo- 
nal ensemble (GOE) of random matrices to a very high 
degree of accuracy. Recent studies on periodic, three- 
dimensional, molecular network-forming liquids have also 
shown that the statistical fluctuations obey the charac- 
tersistics of the GOE 2j| at various temperatures and 
for various densities. Even in amorphous alloys (28| . the 
statistical fluctutations have been reported to obey the 
characteristics of random matrix theory. Another univer- 
sal feature that has been reported in earlier and recent 
studies on periodic three-dimensional amorphous systems 
using various model potentials is the density of states 



function approaching a limit that is independent of the 
explicit functional form of the potential in the amorphous 
regime. The reasons for this universality have been sug- 
gested in [ii,lii,|2o|. 



The concepts of localization and delocalization of vi- 
brational states has also been a focus in these studies, 
with delocalized states being associated with GOE statis- 



tics [17|, l27| and localized states resulting in Poissonian 
statistics tl7|. A popular measure of localization is the 



2articipation ratio [Ti, [3, [B, 0, [III, [H, 0, [H S E 
27[ . Values for participation ratios for three-dimensional 
amorphous clusters 20, 21, 2^[2^, network- forming liq- 
uids 2j| and atomic liquids [27'| indicate that the major- 
ity of the states are delocalized, supporting the apparent 
universality of GOE behavior in the statistical fluctua- 
tions. 

The recent years have also seen an increase in the stud- 
ies of vibrational properties of nanocrystalline materials 
both theoretically [29] and experimentally [30| . The new 
experimental techniques and methods have made it pos- 
sible to explore the systems that are very small in size [sij 
and increased computational power has even made it pos- 
sible to explore such systems in reduced dimensions [32| . 

In the present work we numerically investigate the vi- 
brational spectra of two-dimensional clusters. It is well 
known (e.g. in critical phenomena) that dimension has a 
strong effect on universal behaviors and a detailed study 
in a lower dimension would be a severe test of the robust- 
ness of the observed universalities in three-dimensional 
systems. E.g. as suggested in (33| . in quantum mechan- 
ical systems such as two-dimensional kicked rotors for 
lower dimensions {d < 2), all the electronic states are ex- 
ponentially localized, while for higher dimensions (d > 2) 
there are extended as well as localized states. Also, dis- 
ordered molecular systems display novel behavior in two- 
dimensions, e.g the occurance of the hexatic liquid phase 
that has a long range bond order without translational 
order 3J]. Hence, checking the universality of the vi- 
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FIG. 1: (a). A typical starting configuration of a two- 
dimensional cluster with = 2000. (b). Stable two- 
dimensional cluster obtained using method of Homtopy with 
iV = 2000. 



brational spectra for two-dimensional clusters is a severe 
test. Our interest is to see the effect of dimensionality on 
the states being localized or delocalized with the main 
focus being on statistical fluctuations for the present two- 
dimensional case. 



METHODOLOGY AND DENSITY OF STATES 

The model used for the interactions between the par- 
ticles is the Morse potential |3S] . Over the years various 
studies have demonstrated the usefulness of this poten- 
tial in simulations with references [H, [s^] serving as a 
few examples. The form of the potential energy V of the 
system is, 



V = ^{exp [-2a (n^ 

j>i 



l)]-2exp[-a(r,,-l)]}. (1) 



The potential is thus a function of inter-particle distance 
rij and the (positive) parameter a; the sum is over all 
pairs of particles in the system. The a parameter can be 
tuned to fit a variety of systems ranging from metals such 
as sodium with a — 3.15; van der Waals bonded systems 
such as rare gases with a = 6; to the very short-rariged 
interactions of the Cgo molecule with a = 13.62 37 1. 
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FIG. 2: Color online (a). Density of states {G (uj)) vs. nor- 
malized frequencies of normal modes plotted for various val- 
ues of a. Area under each curve has been rescaled to unity. 
N — 4000 has been used in this calculation, (b). Partici- 
pation ratios (pa) vs. normalized frequencies of the normal 
modes for various values of a. (c). log (pa) vs. normalized 
frequencies of the normal modes for a = 6.0 for A'' = 500 and 
N = 2000 particles. 



The values of a that have been used in the present study 
are 3.5, 6.0, 10.0, 13.0 and 16.0. Using these values of 
a for the interaction potential, we first generate stable 
two-dimensional clusters with N = 4000, where N is the 
number of particles in a cluster. For some calculations, 
clusters with N ~ 500 and N = 2000 particles are used. 

In order to generate stable clusters, we begin by initial- 
izing the particle positions within the limiting distance 
r — N /tt on the x-y plane. Fig. 1 (a) shows an exam- 
ple of a typical starting configuration for N = 2000. In 
order to obtain the potential energy minimum for such a 
configuration, we use the homotopy method of minimiza- 
tion 
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39|. According to this method, in order to find a 
local minimum of a function V ^ we minimize in a series of 
steps the quantity 9V +{l — 9) W, where is a suitably 
chosen simple function (e.g. quadratic) and 9 is varied 
from to 1 in a finite number of steps, typically 20. The 
minimized configuration for one value of 9 serves as the 
initial configuration for the subsequent step. Fig. 1 (b) 
illustrates a two-dimensional stable cluster generated by 
this method. 

At each local minimum generated by this method, a 
hessian matrix is constructed using the position coordi- 
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FIG. 3: Integrated density of states for a single spectrum for 
a = 6.0 with = 4000. Data in the region I has been used 
for unfolding and to analyze the statistical fluctuations 



system as per the notation used by Hudon et al [32l |. 
The peaks shift towards right with increasing a and the 
overall structure varies very slowly for higher values of a. 

This is also notably different than the earlier study 
conducted on a periodic three-dimensional system us- 
ing a similar generalized Morse potential [20|, |23i]. In 
the three-dimensional case, it was possible to generate 
states with variable amounts of disorder and the anal- 
ysis of disordered states showed that the two-peak be- 
havior of DOS changed over to a one-peak as the value 
of a was increased and by the time a reached 16, only 
one peak remained. Further increase in a values lead to 
small changes in the overall DOS shape [l^, [lH. How- 
ever in the present case, potential energy minimization 
results in states that are weakly-disordered rather than 
being totally amorphous and at present, it is not clear 
how one could generate stable clusters in two-dimensions 
that have variable amounts of disorder. However a simi- 
lar feature at both two and three dimensional case occurs 
at higher values of a, where DOS function varies very 
slowly with further increases in a. 



Fig. 2(b) shows theparticipation ratios [151, |17|, 
HI, [il S S, m, 0, Hi calculated using the eigenvec- 
tors corresponding to each eigenvalue for all values of a. 
Mathematically, the participation ratio is defined as. 



nates of particles at that minimum and is diagonalized 
using standard methods to obtain the eigenvalues as well 
as the eigenvectors. The eigenvectors correspond to the 
normal modes and the eigenvalues (A) are related to the 
frequencies (w) of the obtained normal modes according 
to uj = VA. The normal modes obtained in this manner 
are also known as quenched normal modes (QNM). A 
quantity that is of central physical interest is the density 
of states (DOS) function, the histogram of values of w, 
denoted by G{uj). In Fig. 2(a), we plot G (w) for several 
values of a. To improve statistics we combine informa- 
tion from the available distinct quenched configurations 
(atleast 75 for each case) by averaging over all the con- 
figurations. After obtaining the raw averaged DOS func- 
tion, the average frequency is normalized to unity and 
G(a;) is rescaled so that the area under the curve is also 
unity. This process enables us to compare the DOS for 
different values of a. Fig. 2(a) clearly shows the exis- 
tence of a two-peak behavior for all the values of a stud- 
ied which suggests that the obtained local minima have 
weakly-disordered domains with grain boundaries. This 
is further evident from Fig. 1(b). This two-peak nature 
observed in the DOS function is also consistent with the 
recent results obtained by Hudon et al. [13] where bulk 
nanocrystalline materials in two-dimensions are studied 
using Lennard-Jones interaction potential and a similar 
form for G{uj) is observed and it also shows that the 
present system behaves as a "normal" two-dimensional 



(2) 



where is the projection of the eigenvector (labelled by 
a) onto particle i. For extended modes, p is of the or- 
der of unity and does not depend on system size, while 
for localized modes it will scale inversely with the system 
size. Calculations for participation ratios have been done 
using N — 500 and N — 2000 particles. For each value of 
a on this plot, information from various quenched con- 
figurations is combined to improve statistics. For all the 
a values studies, maximum values for participation ra- 
tios stay below or close to 0.6 which might suggest that 
more states are localized in two dimensions as compared 
to three-dimensions. However, closer examination of par- 
ticipation ratios plotted on a semi- log scale in Fig. 2(c) 
with varying N shows that the states are definitely not 
localized over the wide spectrum that has been used for 
studying statistical fluctuations. Since the two systems 
differ in size by a factor of 4, localized modes in the larger 
system should have participation ratios approximately 
1 /4 of those in the smaller system. We find that differ- 
ences on this scale only occur towards the higher values 
of Lo. Hence, even though the participation ratios in the 
middle part of the frequency spectrum that is used in the 
analysis of statistical fluctuations are close to 0.5 or 0.6, 
these still have behavior consistent with extended modes 
of the system. 
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FIG. 4: Color Online (a). Probability density p{s) for nor- 
malized nearest neighbor spacing (s) for various values of a. 
Also shown is the prediction for the GOE. A'' = 4000 has 
been used in this calculation, (b). Variance of the number 
of levels in intervals of length r shown as a function of r for 
various values of a. Also shown is the prediction for the GOE. 
N = 4000 has been used in this calculation. 



FLUCTUATIONS 
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FIG. 5; Color Online (a). Skewness parameter of the distri- 
bution of n(r), the number of levels in intervals of length r 
shown as a function of r for various values of a. Also shown 
is the theoretical prediction for the GOE. A'^ = 4000 has been 
used in this calculation, (b). Excess parameter of the distri- 
bution of n{r), the number of levels in intervals of length r 
shown as a function of r for various values of a. Also shown 
is the theoretical prediction for the GOE. N = 4000 has been 
used in this calculation. 



We now investigate the statistical fluctuation proper- 
ties of the DOS. In the present case, the fluctuation prop- 
erties are computed for A's, the eigenvalues. For a par- 
ticular inherent structure, we denote the elements of the 
obtained spectra by X{i) with i = 1,2,..., 2N. Since 
the present system is two-dimensional, the first three el- 
ements in the spectra will be zero and the remaining 
{2N — 3) positive frequencies are characterized by defin- 
ing a mean local density as well as the fluctuations around 
it. The first step in compu ting the fluctuation properties 



is to unfold the data [2l|, 122, 123!, |2J, |4d, |4l|, 14^ . This 



process enables us to transform the eigenvalues in such 
a way that the average spacing between two successive 
eigenvalues is unity. 

For unfolding the data, we define H{X) to be the num- 
ber of frequencies equal to or less than A as shown in Fig. 
3, and >S'(A) be a smooth function that passes through the 
staircase function H in the best-fit sense. In the present 
case there is no single function that passes smoothly 
through the whole of H in the best-fit sense. We leave 
out small regions (about 10% of the levels) from each 



end of the spectrum and for the remaining levels in the 
central region marked as I in Fig. 3, we use a quadratic 



polynomial D{X) 



6 A -|- c A as an approximation 



for S{X) over the remaining spectral region. It must be 
stated that there is no particular reason to leave out lower 
end region of I other than obtaining better fitting to the 
data however, it is evident from Fig. 1(b) and Fig. 1(c) 
that towards the upper end of I, more states are local- 
ized and hence for the analysis of this paper, it is the 
central regime that we are interested in. The values of a, 
b, and c are obtained by a standard least-square fitting 
procedure. We also calculate the misfit function 22, 23[ 
corresponding to each of the fits to check how well D{X) 
approximates S{X). The plot of the misfit function we 
find is qualitatively similar to the plot in Fig. 1(b) of [2^ . 

In order to reduce the mismatch to lower degree, we 
eliminate subregions where the mismatch function has a 
very irregular behavior. In the remaining regular subre- 
gions we fit quadratic function to the misfit function, and 
we correct D{X) by these quadratic functions to obtain 
the desired unfolding functions. For each of the fluctu- 
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ation properties reported in this paper, we combine the 
data from aU the subregions of all the spectra. 

Note that this process of unfolding is quite different 



from the one used in references [211, [22, |23( in which an 
exponential function has been to unfold the data followed 
by the quadratic correction. In the present case, we could 
not find any suitable S{X) that could pass through H. 
Hence we have to resort to this method. However, the 
advantage of this method of unfolding is that it provides 
a suitable way of handling with more complex functions 



24| 



that have no well-defined analytical forms 

The first fluctuation property that we report here is 
p(s), the distribution of the normalized nearest-neighbor 
spacings s of the frequencies of the unfolded spectra. For 
this, we first complete the process of unfolding described 
above for each of the spectra using D{X). At this level, 
we apply the quadratic correction to this D(X) ignoring 
small regions where the misfit function is very irregular. 
In this way we obtain the unfolded spectra for each of the 
spectra. Selection of a random individual spectrum and 
analysis of each of them separately indicates that the the 
spectrum has fluctuation properties associated with the 
Gaussian orthogonal ensemble (GOE) of random matri- 
ces. To improve statistics, we combine the data for all 
the four regions of all the unfolded spectra. This has 
been plotted in Fig. 4(a) for three values of a along with 



the theoretical prediction [40, |41|, [42, |43( . It can be seen 
that the agreement with the theoretical prediction is ex- 
tremely close, in spite of the values of the participation 
ratios being lower than those found in three dimensions. 

The second fluctuation property we report is the quan- 
tity (r), the variance of the number of levels n (r) 
within an interval of length r located randomly in the 
unfolded spectrum. This is plotted in Fig. 4(b). It 
must be emphasized that the quadratic correction ap- 
plied to the fitting function D (w) is very important in 
calculating (r). This calculation is extremely sensi- 
tive even to very small errors in the approximation to 
S. The contribution of any such error to (r) grows as 

, whereas the GOE prediction for (r) grows only as 



In(r') [21|,|22,|23|. Values for a = 3.5 almost overlap to the 
theoretical prediction whereas we do see a shift for the 
other two a values. This may be due to the following pos- 
sible effects: (i) As explained in 2l|, [22, l23[ , the exact lo- 
cations of the irregular regions vary in the contour of the 
misfit function and this might add to the observed shift. 
A detailed analysis using a spectrum — specific [i^, [23l | 
choice of subdomains might help in determining the 
strength of this effect on the observed deviation, (ii) As 
observed from Fig. 1(b), values of participation ratios 
decrease with increasing a and this decrease might have 
a role to play in the observed deviation. 

In Figs. 5(a) and 5(b), we plot the skewness and ex- 
cess parameters [13] of the fluctuations. Also included 
in these plots are the predictions for the GOE. These 
predictions have been calculated on the basis of a large 



ensemble of 500 x 500 matrices belonging to the Gaussian 
Orthogonal Ensemble. Again we observe that the agree- 
ment with the theoretical prediction is extremely close. 
Plots show only one value of a just for clarity. Other- 
values of a also have the same level of agreement to the 
theory as a = 6. 



CONCLUSIONS 

Our results for the DOS of two-dimensional clusters 
shows the survival of a two-peak behavior when the width 
of the potential well describing the Morse potential is re- 
duced progressively thereby indicating a weak-disorder in 
contrast to the three-dimensional case. The participation 
ratios suggest that the vibrational spectrum has behavior 
consistent with the extended modes of the system. Fur- 
ther, the vibrational spectra has fluctuation properties 
associated with the GOE of random matrices. For each 
of the fluctuation properties, agreement with the GOE 
prediction is extremely close. The observed shifts in the 
case of the E^ (r) calculations suggest that magnitude of 
the participation ratio might have a role to play in this 
behavior, however, since calculations of the normalized 
nearest neighbor spacings distribution and skewness and 
excess are extremely close to the theoretical prediction, it 
may not be incorrect to say that the system still follows 
GOE. 

The work leaves us with the following challenges as 
a part of the future work: (1) The first challenge is to 
generate local minima that have a variable degree of dis- 
order in two-dimensions, at this moment it is not clear 
how to generate such states in two-dimensions that would 
allow us to do a systematic study with respect to disor- 
der. (2) What should be the level of disorder before we 
see a convincing departure in the statistics from GOE. 
The present work shows that the statistics are of GOE 
type even in a weak-disorder. The reasons for this effect 
are yet to be determined. (3) As mentioned by Hudon 
et al. [3^ in their studies on bulk nanocrystalline ma- 
terials that density can influence thermal properties of 
nanocrystalline materials. Using this argument, it would 
of great interest to generate local minima as a function 
of densities and repeat the calculations for statistics. 
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